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Abstract. 

Resonant solutions of the quantum Schrodinger equation occur at complex 
energies where the S-matrix becomes singular. Knowledge of such resonances is 
important in the study of the underlying physical system. Often the Schrodinger 
equation is dependent on some parameter and one is interested in following the 
path of the resonances in the complex energy plane as the parameter changes. 
This is particularly true in coupled channel systems where the resonant behavior 
is highly dependent on the strength of the channel coupling, the energy separation 
of the channels and other factors. 

In previous work it was shown that numerical continuation, a technique 
familiar in the study of dynamical systems, can be brought to bear on the problem 
of following the resonance path in one dimensional problems [lj and multi-channel 
problems without energy separation between the channels [2]. A rcgularization 
can be defined that eliminates coalescing poles and zeros that appear in the S- 
matrix at the origin due to symmetries. Following the zeros of this regularized 
function then traces the resonance path. 

In this work we show that this approach can be extended to channels with 
energy separation, albeit limited to two channels. The issue here is that the 
energy separation introduces branch cuts in the complex energy domain that 
need to be eliminated with a so-called uniformization. We demonstrate that the 
resulting approach is suitable for investigating resonances in two-channel systems 
and provide an extensive example. 
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1. Introduction 

The presence of resonances can drastically increase the yield of many quantum 
mechanical reactions [3]. A resonance is an intermediate state that is formed when 
reagents collide with an appropriate energy and form an excited complex that decays 
into reactants or products. In molecular reactions these intermediate states are often 
electronic excited states. In the Born-Oppenheimer picture the molecular dynamics 
temporarily follows the electronic potential energy surface formed by the resonant 
state. 

This picture forms the basis of state-of-the-art ab initio calculations where, first, 
the position and lifetime of the resonance is determined using electronic scattering 
calculations with the nuclei fixed in space. This calculation needs to be repeated for 
every possible configuration of the nuclear positions probed by the chemical reactions 
and this results in the resonant potential energy surface. In the second step, the 
nuclear dynamics is simulated on the resonant potential energy surface that leads to 
the reactants. This approach has been successfully used to calculate the yields of 
processes such as dissociative electron attachment to the water molecule [4j [5j [6] and 
vibrational excitation of carbon dioxide [7J [8] , to name a few of the processes that are 
mediated by a resonant state. 

This work focuses on the first step outlined above, where the potential curves are 
calculated as a function of the nuclear degrees of freedom. We will study a model 
consisting of a coupled Schrodingcr equation and construct, in an automatic way, the 
potential curve of a resonance that becomes bound as the parameters in the equation 
arc varied. The method tracks the resonances accurately, even in parameter ranges 
where the resonance is too broad to be tracked with traditional methods such as 
complex scaling [9]. 

In this article, as is regularly done in the literature, a resonance energy is defined 
as the complex energy at which the 5-matrix has a pole [TQl [TTj . The real part of the 
pole position defines the physical resonance energy and the imaginary part defines the 
resonance decay width (or inverse lifetime). The advantage of this approach is that 
a bound state energy is also a pole of the S'-matrix, albeit with zero decay width i.e. 
infinite lifetime. As such, bound states and resonances are both treated in the same 
way which makes it more convenient to trace resonant states as they become bound, 
or vice- versa, when a problem parameter changes. The poles of the S'-matrix will be 
found by a Newton iteration applied to a function that is proportional to the Jost 
function. 

We aim to trace these states with numerical continuation. It is a technique that 
has found widespread application especially in the dynamical systems community. 
Assume one is given a solution (uq, Ao) of a set of n nonlinear equations F(u, A) = 0, 
where F : R™ +1 — > R™. A second solution (ui,Ai) is then constructed numerically 
by applying a predictor-corrector scheme. Repeated application constructs an 
approximation of the implicitly defined solution set of F. One of the well-known 
numerical continuation techniques is pseudo-arclength continuation proposed by 
H. B. Keller [12]. It has been implemented in computer programs e.g. AUTO |T3l [14] . 
LOCA |15] and other numerical continuation libraries. Because the corrector step in 
these methods is based on Newton iterations, the derivatives of the function F should 
be Lipschitz continuous to guarantee fast convergence. 

A traditional method to find the resonance position and width is complex 
scaling [9]. In this method, one applies a complex scaling transformation, r — > r 10 , 
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to the reaction coordinate. This turns the resonance wave function into a square 
intcgrable function. After the transformation the resonant state is then part of the 
discrete spectrum of the Hamiltonian. The method, and its exterior variant (ECS) 
which has the advantage that it leaves the interactions in the inner region unchanged, 
have been successfully applied to find resonances in molecular systems such as HCO 
[16] . NelCL |17| and in many other examples in atomic, molecular and nuclear physics. 

Complex scaling, however, has its limitations. First, in numerical calculations 
the resonance position depends slightly on the choice of the rotation angle 8. This 
is documented for example in [5]. A second limitation is that only resonant states in 
a limited region in the complex energy plane can be found, in particular, a pie slice 
between the continuum spectrum, which is rotated 29 downwards from the real axis, 
and the real axis. Virtual states for example, have a purely imaginary wave number 
and lie outside this region. As such they are hard to find with this method. 

A resonance often transforms, very shortly, into a virtual state before it becomes 
a bound state as the parameters of the system change. To understand how resonances 
and bound states are connected through these states it is necessary to have a 
mathematical description that can handle these virtual states. 

This paper investigates the application of numerical continuation to trace resonant 
states in a coupled channel Schrodinger equation as the parameters of the problem 
change. Unfortunately, the application of numerical continuation to trace resonant 
states is not without challenges. For Newton's method to work efficiently it requires a 
function whose derivatives are Lipschitz continuous |18| . The S-matrix does not satisfy 
these smoothness conditions. In particular when a resonance makes the transition to 
a bound state, a pole and a zero of the S-matrix meet and straightforward application 
of numerical continuation fails to trace the zero of 1/S. 

In [1] it was found, for one-dimensional quantum systems, that it is possible to 
apply the continuation to a regularized function derived from the S'-matrix because 
that function does satisfy the necessary smoothness conditions. For several realistic 
potentials the resonances were tracked as parameters in the system were varied. 

In [5] the method has been extended to one-dimensional, many-channel problems 
where all channels have the same asymptotic energy threshold. Applications have 
demonstrated the viability of the approach in tracing the parameter dependence of 
the resonance in the system. However, as it stands, the method does not apply to 
systems where the channels have unequal energy thresholds. In this case, branch cuts 
occur in the complex energy plane making the method invalid. 

This is a significant restriction for certain application areas, e.g. most problems 
in molecular dynamics have such unequal energy thresholds. In this paper, we will 
investigate this difficulty and show that in the two-channel case it can be overcome 
by introducing an appropriate uniformization of the complex plane. In addition we 
demonstrate the numerical practice of the method using an extensive example. 

We proceed as follows: In section [2] we formulate the equations of interest 
in the context of quantum mechanical systems. Section [3] describes the so-called 
scattering matrix or S'-matrix, its relation to resonances and bound states as well 
as the complex geometries that occur in the S-matrix of coupled channel systems. 
Section |4] gives a brief overview of the background and applications of numerical 
continuation techniques. Section [5] outlines the details of our implementation of the 
methods described in the preceding sections. In section[S]we present an excerpt of the 
results we have obtained and compare with the ECS method and in section [7] we give 
an outlook for possible applications and future studies. 
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2. Coupled channel Schrodinger equation 



The time independent Schrodinger equation for TV coupled channels with a spherically 
symmetric potential reads 



1 d 2 



-I 



L(L + I) 



V(r, A) + S *(r; E, A) = E*(r; E, A), (1) 



2/i dr 2 2fir 2 

where (r £ M + ) is the radial coordinate, I is the N x N identity matrix, S is the 
diagonal matrix of channel thresholds L is the diagonal matrix of channel angular 
momenta li, V(r, A) is the matrix of channel and coupling potentials Vij(r, A) which 
depend on a problem parameter A € K and, finally, \t(r) is the matrix of channel wave 
functions. Depending on the properties of the potential matrix V(r, A) at infinity the 
behavior of the solutions \&(r) will differ significantly. We assume so called short range 

oo and is less singular than r~ 2 

and an incoming 



as r 



interactions: Vij(r) vanishes faster than r 
in the origin r = 0, see also 

Given homogeneous Dirichlet boundary conditions at r 
plane wave, the asymptotic (r — > oo) solutions behave as: 



*(r;E,X) = - 



h£(Kr) - h+(Kr)K-zS(E, \)K* 



(2) 



where K 



y // 2fi(EI — S) is the matrix of channel momenta ki 



h±(Kr) 



v / 2/^(^ - ^) and 



(3) 



V 



h l N ( k Nr) 

is the matrix of spherical Riccati-Hankel functions associated with the various channels. 
The first term in ^ is the partial wave expansion of the incoming plane wave, the 
second term represents, for each incoming partial wave, the outgoing wave in all the 
channels. The matrix S(E, A) is the so-called scattering or S'-matrix. 

Our main object of interest is S(E, A) because it contains all the information about 
the scattering process. It can be obtained from the coupled channel wave functions \l/ 
through the expression: 

S{E,X) = K-'W \hl(Kr ),*(r ;E,X)\ W \h+(Kr ), *(r ; E, A)l K 
where the W stands for the Wronskian of two functions, whose usual definition 



W[f(x),g(x)]=f(x 



dg(x) df(x) 



dx 



dx 



is extended to matrices of functions as in [19]: 

dB(x) 



W[A(x),fl(x)]=^ (x) 



dA T (x) 



B(x), 



(5) 



(6) 



dx dx 

where T is the transpose of a matrix, which vanishes in the specific case of 
expression (@|. 

The Wronskians in equation ((4]) are evaluated at a point ro outside the range 
of the potentials, usually near the edge of the computational domain. Therefore, 
given a numerical method for solving the Schrodinger equation (TTJ), the S'-matrix can 
be obtained numerically by evaluating matrix expression Q which also involves the 
derivative of the wave functions. The details of our numerical implementation are 
given in section [5j 



Numerical Continuation of Bound & Resonant States 



5 



3. Extracting the resonant and bound states from the S-matrix 

We are interested in solutions of equation ([T]) that correspond to the bound states and 
the resonances of the system. Many characterizations of these special solutions exist, 
yet the theoretically fundamental definition describes both bound and resonant states 
as an eigenstate of equation ([T]) with purely outgoing wave functions at infinity as the 
second boundary condition [TU1 [11] . 

A consistent alternative definition for multichannel systems interprets these states 
as having energies E for which det(S(£ l , A)) exhibits a pole. This is consistent with the 
first definition, as is easy to see when one looks at the expression for the asymptotic 
wave function ||2J). If dct(S(_B, A)) has a pole, then in at least one eigenchannel (i.e. the 
channels defined by linear transformation to diagonalizc the S(E, A) matrix) a resonant 
solution proportional to a purely outgoing wave occurs. A more formal approach relies 
on the introduction of Jost functions, for which we refer to the literature [10| ITT] . 

From the definition of the bound and the resonant states as a pole of dct(S(£', A)) 
it is now possible to study the evolution of these states in terms of a changing system 
parameter A £ 1. As A traverses the parameter space, the bound and resonant states 
change position and lifetime. It is of interest for many applications to know the explicit 
dependence on the parameter of choice because it leads to the resonant and bound 
state potential surface. 

It would be straightforward to define resonance trajectories as curves in the 
complex _E-plane parameterized by A i.e. {(E, A) G C x K | dct(S(i?, A)) -1 = 0}. 
Although this is theoretically correct, it is not feasible numerically because of two 
reasons. First, in the general case, the coupled channel S-matrix is a multi-valued 
complex function of the energy and "lives" on a multi-sheeted Riemann surface with 
branch cuts for every threshold value. Second, near threshold parameter values where 
bound and resonant states meet, multiple poles and zeros of the S'-matrix coalesce, 
thereby destroying local smoothness properties. This is a consequence of well known 
symmetry properties of the S'-matrix [TjJJ [TTJ . 

The presence of branch cuts and the coalescence of zeros and poles have a 
strong negative impact on the convergence of the Newton iteration used in solving 
dct(S(i?, A)) -1 = when calculating the resonance trajectory. A study of specific 
cases will give us insight in how these issues can be addressed. 



3.1. The case of N channels, equal thresholds 



This case has been investigated in [2] and we briefly review the results. In eq. (HJ the 
matrix of channel thresholds H is zero (i.e. we take the threshold energy to be zero 
in all channels). Consequently, the matrix of channel momenta K becomes kl with 
k = \J2[iE. This simplifies expression ([2]) to 



h-(fer) - h+(kr)S(E,X) 



such that dU) becomes 

S(E,X) = w[h£(fcro),*(r ;-B,A) 



w 



h+(kr ),*(r ;E,\) 



(7) 



(8) 



Because k = ^2/iE, the function S(E, A) is defined on a two-sheet Riemann surface 
with a branch cut on the negative real axis E £ Mr . This destroys continuity of the S'- 
matrix near the branch cut and makes numerical continuation difficult. Fortunately, 
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in this case, one can express the jfj-matrix in terms of fc by introducing S(fc, A), a 
single-valued function of k. Afterwards, the results can be easily translated back 
to the complex iJ-plane with E = |-. In mathematical terms such a procedure of 
rcparamctrization of a multi-valued complex function to a single-valued function is 
referred to as a uniformization. 

Another issue that has to be dealt with is the coalescence of poles and zeros 
near the k = threshold value. This can be circumvented by effecting the numerical 
continuation on the function 

F(k,X) =dct (/C(S(fc,A) -I)" 1 ) , (£)« = k 2U+1 (9) 
= (j[k 2k+1 ) /det(S(fc,A)-I), (10) 

instead of directly on the 5-matrix. In (2] it is shown that this procedure does not 
introduce false solutions, i.e. the zeros of F are precisely the poles of S, and that it 
eliminates any singularities at k = 0. This is an extension of a similar result for the 
one-dimensional, single-channel systems case [T] . 



3.2. The case of two channels, unequal thresholds 



To highlight the difficulties that arise in the case of unequal thresholds we focus on a 
two channel model. Each threshold £j introduces a separate branch cut (— oo,£j], see 
figure 1(a) Rcparamctrization of S in terms of either of the channel momenta hi does 
not provide a viable uniformization. In each of the ki or fc2-pla-n.es there is a branch 
cut at [—b, b] where b = y/2fi(£ 2 ~ £1) (without loss of generality we take £1 < £2), as 
illustrated in figure 1(b) The cut disappears only when the thresholds coincide. 

For two channel systems a uniformization exists, given in |10j and modified slightly 
in [19]. One introduces u £ C such that 



E(u) 



1 



2u 2 



The corresponding expressions for channel momenta fci and fc2 are: 

^7« 2 -i 



fci(u) 



k 2 (u) 



6 



/'- 



6 - a u' + 1 



/'- 



2 u ' 

and the 2x2 S'-matrix is now explicitly written as 



where 



S(u,A) 



W±(u,A) 



fci(u) 
k 2 {u) 



h h (ki(u)r ) 




W_(u,A) [W+^A)]- 1 



(11) 

(12) 
(13) 



ki(u) 
fc 2 (u) 



,(14) 



hf(k 2 (u)r 




±i>u{r ;E(u),\) f^ l2 {r Q ;E{u),X) 
frhi(r ;E(u),\) ±i> 22 {ro;E(u),X) 











if>n{r ;E(u),\) ipi 2 (r ;E(u),X) 
il>2i(ro',E(u),X) ^ 22 (r Q ;E(u),X) 



where Vij( r o; E(u), A) stands for the (z,j)-th element of the matrix wave function \& 
calculated for a specific parameter A and evaluated at a point tq. 
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S*(*a) 






-6 









" 5R(fc 2 ) 



(a) With unequal thresholds, (b) Even in the ki or fc2-plane 
branch cuts at (— oo, there is a branch cut at [—6,6] 



Figure 1. Branch cuts of S in E- and fc-representation for a two channel system. 
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Figure 2. Uniformization of the four-sheeted Riemann surface of S(-E) using 
the transformation Hilt unravels the sheets onto four different regions of the 
complex u-plane. The blue region (+, +) corresponds to the physical sheet 
|u| > 1,3?(m) > 0. Bound states (BS) E G (— co,£i] lie on u G [l,+oo). Scattering 
states (Cont.l) E g [£1,62] map to 11 = e iB ,0 e [— f ,0]. Scattering states (Cont.2) 
E 6 [£2, +00) map to u S [— i, —ioo). 

The effect of this uniformization can be visualized by looking at the different 
parts of the u-plane and how they map on the i?-plane as illustrated in figure [2] 
Four regions (±, ±) in the u-plane can be distinguished. They are separated by the 
imaginary axis and the unit circle and correspond to the four different sheets of S(E, A) 
in the Emplane. The region labeled as (+, +) bounded by |u| > 1 and 5i(u) > is 
mapped to the physical sheet. Physical bound states E € (—00,^1] are located on 
u e [l,+oo). The physical continuum E £ [£1,^2] ma P s to the quarter unit circle in 
the fourth quadrant u = e l9 ,9 £ [— f, 0]; whereas the physical continuum E £ [£2, +00) 
is mapped to u £ [— i, —ioo). A more detailed description of the different regions in 
the u-plane can be found in (T9l [10] . 

Continuation paths can traverse different regions in the u-plane, and as such, 
different sheets in the complex energy plane. The focus of this work is concentrated 
on finding these trajectories, independently of their precise physical meaning. The 



Numerical Continuation of Bound & Resonant States 



8 



results wc obtain are presented in the u-plane. However, they can be translated to 
the complex energy plane if one is interested in physically observable quantities. 

Note that this uniformization is strictly limited to a two channel case. A similar 
procedure for three channel systems is much more involved, as indicated in [10] , To 
the knowledge of the authors there is no generalization for N channels. For a thorough 
discussion of the analytical properties of the S'-matrix and related functions in many- 
channel problems we refer to |10| 111] . 

Having established a feasible way of extracting the S-matrix for a specific u E C, 
we are still faced with the problem of coalescing poles and zeros as discussed in the 
previous case of equal thresholds. Fortunately a straightforward extension of the 
regularization procedure (|10[) can be applied. Similarly, we create the function 



F( w ,A) = det(/C(S( M ,A)-I)- 1 ), {JC) U = k* h+1 (u) (17) 
= (f[kf I+1 (u)) /dct(S( u ,A)-I), (18) 



\i=l 



which in our case reduces to 

F{U ' X) - dct(S( u ,A)-I 2 ) ■ (19) 

This is the function whose solution set will be approximated through numerical 
continuation to obtain the trajectories of the resonant S'-matrix poles. 



4. Numerical continuation of resonances 



In the previous section we have identified a resonance trajectory starting at (E ,Xq) 
as the implicitly defined curve E = E(X) with 

F(E(X),X)=0, E(X )=E , (20) 

where the function F is defined in equation (|19j). We are now interested in constructing 
such trajectories automatically and robustly. 

In dynamical systems similar equations arise in the study of steady states of 
parameterized ODEs and efficient methods have been developed to find the solution 
curves in terms of varying parameter values. In this context, one is interested in 
finding the solution of an under determined system of nonlinear equations 

F : — >R n : x = (u, A) i — ► F(x), (21) 

connected to an initial point xo = (uo, Ao). 

Many of these problems are computationally intensive and efficiency is a key 
concern in the numerical studies. In particular, the number of evaluations of the 
function F should be kept to a minimum. In addition, the solution components often 
have complex geometries with intersections and bifurcations. The study of bifurcations 
generally involves rigorous stability analysis of the underlying solutions and is a 
complicated subject on its own. In this work only the so called "simple" bifurcation 
points can occur. These manifest themselves as two intersecting solution branches and 
are characterized by the dimension of the null-space of the Jacobian F x (x t ) being 2 
in a point x t , which is called a "branch point". In many problems, however, the 
bifurcations are much more involved and a thorough treatise on bifurcations and 
stability of solutions can be found in [12j [21] [22j [23] . 
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Figure 3. Schematic representation of pscudo-arclcngth continuation. 

Numerical continuation is the process of solving the equation pip by constructing 
successive approximate solutions on the path starting at the known solution Xrj. A 
good overview of these techniques can be found in [23] . There are essentially two 
different approaches to tracing such paths viz. piece- wise linear methods and predictor- 
corrector methods. We will use one of the latter methods. Typically, one first makes 
a predictor step (Euler prediction is commonly used) that estimates the next point 
by following the tangent tj to the curve at the current point Xj for a certain small 
distance As as in x^ +1 = Xj + Ast,. Next, a corrector step is applied to converge to 
a solution x^ +1 — > x.; + i. Quite often Newton iterations are used as a corrector. 

One such predictor-corrector method, the one we will use in this paper, is the 
pseudo-arclength continuation. Its corrector step consists of Newton iterations on the 
system -F(x) = augmented with an additional equation that constrains the iterations 
to a hyper plane through x^ +1 and perpendicular to the tangent thereby giving the 
next point on the solution curve (see figure 

A robust implementation of numerical continuation that can detect branch points 
and continue branching solution curves is provided by AUTO [13[ 114] , which we have 
used in this work. Other well-known implementations of numerical continuation 
algorithms include LOCA [15] which is a part of the Trilinos framework [24], 
Matcont [35] (a Matlab implementation) and Multifario [53] which allows multi- 
parameter continuation. 

5. Numerical implementation 

Due to its ability to automatically find and follow bifurcating branches, we have opted 
to use AUTO for numerical continuation. Therefore the main component that has to 
be provided to implement methods described in previous sections is the actual code 
that calculates the numerical values of the function F(u, A), see eq. (fT9]) . for a given 
« 6 C = I 2 and A G K. 

We have solved equation |T]) both with a 0(h 5 ) renormalized Numerov method 
detailed in [37] [25] ■ We also need the derivative of the numerical wave function at the 
end of the computational domain for a specified energy E{u) in order to calculate the 
Wronskians in (fT4]) . We compute this derivative using a technique also described in 
[28] and which is of order 0(/i 4 ). Finally, we need to calculate the function F given 
by expression (fl"9]) . 

We emphasize that the numerical continuation results are independent of the 
underlying numerical solver of the Schrodinger equation and that both 0(h 5 ) and 
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the 0(h 2 ) methods give the same qualitative results. However, we have found that 
higher order methods for both the wave function and its derivative lead to significantly 
more robust continuation curves and allow for wider energy ranges. For this reason 
we selected the higher order renormalized Numerov method to generate results of 
section [6j 

This approach has been implemented in CH — h and is used as a driver routine by 
the AUTO program. 

6. Examples and results 

As an example of the complicated geometries this method is able to cope with, we 
consider a coupled channel s-wave {l\ = \<x = 0) system with Gaussian potential wells 
both as the channel potentials and the coupling. The 2x2 potential matrix has the 
elements: 

Vu(r,Xi) =-\ ie -- ^ 
Vi^j(r, A c ) = \ c er r , 

where Ai, A2 and A c denote the potential strength of the first, second and coupling 
channels respectively. The channel thresholds are chosen £1=0 and £ 2 = h and the 
mass is fj, = 1. Equation (TTJ) was solved using the previously mentioned renormalized 
Numerov method on the domain r £ [0,4.8] with 4096 grid points. All of Ai, A2 and 
A c will be used as variable system parameters. We will indicate clearly which of those 
are fixed and which are used as continuation parameters. 

In the uncoupled case (A c = 0), setting A; = 4 gives a system with two bound 
states in each channel whose values can be found in the upper part of table [TJ Using 
these values as starting points we carry out the continuation in terms of increasing 
channel coupling A c while keeping Ai = 4 fixed. The results in the u-plane are 
shown in figure [4] Points corresponding to coupling values 0, 0.2, 0.3 and 0.5 are 
highlighted in the figure. Their corresponding numerical values in the u- and S-plancs 
are summarized in table [TJ Notice that for A c = 0, every value of E is associated with 
two different points in the it-plane which are located on different sheets. Following 
the paths of such two values eventually gives rise to different values in the _E-planc. 
Therefore, one needs to keep track of all of them. We give appropriate labels, see 
table [TJ to distinguish different points in the it-plane. 

While continuing in A c , only slight variations in state's energy occur, yet the 
seemingly minor coupling has profound effects on the evolution of these states in 
terms of potential strengths Ai and A2 while keeping the coupling strength constant. 

To illustrate this behavior we fix the coupling strength and use the corresponding 
u-values of the states as starting points for a continuation in terms of variations of 
both potential strengths Ai and A2 simultaneously. As \ decreases, we expect these 
bound states to move into the resonant regime and influence each other due to the 
coupling. 

Although we have performed numerical continuation starting from all states for 
all subsequent coupling strengths from table [TJ we focus on the case A c = 0.5 as a 
highlighted example. 

The resulting continuation curves are shown in figure [5] as projections on the 
planes 3?(it) x Aj, x A; and 3?(u) x Complex geometries and intersections are 

constructed automatically by the continuation method. An attempt to present those 
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State 


label 








S(ii) 




»(E) 




9(E) 




ch 1, ng 


cin 


-2 


2983975e-01 







-2 


1228484e+00 











cin 


4 


3508575e+00 







-2 


1228484e+00 









m 


cifii 


-4 


5199837e-01 







-3 


8737558e-01 







U 




cini 


2 


2123974e+00 







-3 


8737558e-01 







ch 2, no 




2 


5892712e-01 







-1 


6228484e+00 











C2no 


3 


8620906e+00 







-1 


6228484e+00 









ni 


C2fii 


8 


8019950e-01 


4 


7460388e-01 


1 


1262442e-01 











C2ni 


8 


8019950e-01 


-4 


7460388e-01 


1 


1262442e-01 







A c 


State 


label 




U{u) 




9p(«) 




Sft(E) 




3(E) 




ch 1, no 


cin 


-2 


2923691e-01 







-2 


1352756e+00 











cin 


4 


3623083e+00 







-2 


1352854e+00 









ni 


cifii 


-4 


5179967e-01 







-3 


8789140e-01 







. 2 




cini 


2 


2141945e+00 







-3 


8832855e-01 







ch 2, no 


C2U0 


2 


5963744e-01 







-1 


6127068e+00 











C2"-0 


3 


8517883e+00 







-1 


6129594e+00 









ni 


C2ni 


8 


7757633e-01 


4 


7363933e-01 


1 


1278823e-01 


1 


1579542e-03 






C2ni 


8 


7757633e-01 


-4 


7363933e-01 


1 


1278823e-01 


-1 


1579542e-03 


Ac 


State 


label 








3(u) 




5R(E) 




3(E) 




ch 1, no 


cin 


-2 


2852083e-01 







-2 


1501654e+00 











cino 


4 


3759865e+00 







-2 


1501849e+00 









ni 


cifii 


-4 


5155161e-01 







-3 


8853641e-01 











cini 


2 


2164315e+00 







-3 


8951600e-01 







ch 2, no 


C2no 


2 


6048642e-01 







-1 


6006947e+00 











C2«0 


3 


8395472e+00 







-1 


6012444e+00 









ni 


C2fii 


8 


7428785e-01 


4 


7241720e-01 


1 


1298423e-01 


2 


6183712e-03 






c 2 ni 


8 


7428785e-01 


-4 


7241720e-01 


1 


1298423e-01 


-2 


6183712e-03 


Ac 


State 


label 




K(u) 




3(u) 




5R(E) 




3(E) 




ch 1, no 


cifio 


-2 


2645171e-01 







-2 


1939897e+00 











cin 


4 


4159879e+00 







-2 


1940286e+00 









ni 


cifii 


-4 


5076010e-01 







-3 


9060199e-01 







0.5 




cini 


2 


2235200e+00 







-3 


9328811e-01 







ch 2, no 


C2"0 


2 


6297322e-01 







-1 


5661803e+00 











C2no 


3 


8041557e+00 







-1 


5675876e+00 









ni 


C2fii 


8 


63688796-01 


4 


6837873e-01 


1 


1354314e-01 


7 


3933316e-03 






C2ni 


8 


6368879e-01 


-4 


6837873e-01 


1 


1354314e-01 


-7 


3933316e-03 



Table 1. Numerical values of the states in the u- and E-planes. Different states 
with coupling values (top to bottom) A c = 0, A c = 0.2, A c = 0.3 and A c = 0.5 are 
shown. The colored labels are used in figures to make a clear distinction between 
states. 



schematically is shown in figure[7]as a connectivity graph. In the range Aj € [0, 4] four 
branch points can be distinguished with values detailed in tabled] A three-dimensional 
overview of the continuation curves is shown in figure [6] 

We have carried out the same procedure for four different values of the coupling 
strength using starting values summarized in table [TJ As the coupling between the 
two channels increases, resonant trajectories undergo major qualitative changes. See 
figure IHl for a short comparison of the continuation paths projected on 5ft(u) x Aj. 
The interpretation of those is beyond the scope of this paper although an important 
effect can be observed. In figure [10] a close-up view of two resonant trajectories is 
shown. In the uncoupled case (figure 10(a) ) two independent trajectories of ciUq (and 
C2^o) and Cini are shown. Although the green and red curves intersect they do not 
share a branch point. The addition of a nonzero coupling in figure [lO (b) | changes this 
situation drastically: two more bifurcation points appear and the red and green curves 
are now fully connected (dashed green/red lines). As the coupling strength increases 
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0.5r 

0.4- 
0.3 - 



0.4 
0.3 



Re(u) 

(a) 5R(«) X A c projection 



lm(u) 

(b) £j(u) X A c projec- 
tion 



2 
Rem) 



(c) SR(n) X S(u) projection 



Figure 4. (color online) Continuation paths of bound states of the example 
system in terms of increasing channel coupling strength A c projected on three 
different planes. As the coupling strength increases, slight repulsion of the states 
in both channels can be observed. Circled values indicate starting points used 
for continuation in channel strengths \{ . The values of these points are given in 
table [T] In figure [4(c) | the dashed line represents the unit circle. 



in figure 10(c) the two additional branch points collide and disappear. For an even 
higher coupling (figurc [lO(d)[ ) one can clearly see how the curves have rearranged their 
connections: The cini state is now connected with C2ho through the dashed green/red 
line, whereas c^no shows a connection with the end point of cini from the uncoupled 
case. 

This short example highlights the ability of continuation methods to deal with 
subtle and complex connections. 



6.1. Comparison with Exterior Complex Scaling 

To highlight the advantages of the numerical continuation method applied to the 
function F(u,X), we compare its results with those of a calculation with exterior 
complex scaling (ECS) [9]. This is done for various choices of the parameter A^. First, 
we translate the curves that were obtained by numerical continuation in the u-planc 
back to the i?-plane. These are shown as the green and blue curves in figure [8] For 
clarity, the colors are identical to those in the figures depicting the u-plane. Note 
that two different blue curves map to the same region in the .E-plane. The vertical 
axis is the strength A^ of the channel potentials, while the bottom plane shows the 
complex energy of the resonant state. The real part is the proper resonance energy and 
the imaginary part is the inverse lifetime of the resonance. We have found both the 
exponentially decaying and exponentially growing states with a negative, respectively, 
positive imaginary energy. 

A first bound state is formed as A^ increases and the potential becomes stronger. 
It starts out as a virtual state on the negative real energy axis and, with increasing 
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Figure 5. (color online) Different projections of the continuation curves (A c = 
0.5). The four branch points are highlighted and labeled accordingly. Their 
numerical values are summarized in table [2] The connectivity graph in figure [7] 
gives a clearer overview of the associated connections. 



A,;, its real part increases, goes through zero and becomes negative again. This curve 
is shown in blue in figured] 

For potential strength A; between 0.33538 and 1.54368 we also have a resonant 
state with the real part of the energy between the values of the two thresholds. This 
state is a Feshbach resonance related to the second channel that decays through 
coupling with the first channel, which is open. This state also originates as a virtual 
state on the negative real axis at small potential strengths, similar to the state 
discussed above. However, it is now shifted up by 0.5, the threshold of the second 
channel. Furthermore, because of the coupling to the open channel this virtual state 
has a finite lifetime as soon as it lies above the threshold of the open channel. When 
we increase the potential strength above 1.55204 this resonance becomes a second 
bound state, but before it becomes bound, there is a bifurcation with a virtual state 
at parameter strength A,; = 1.54368. At this bifurcation point the state has negative 
real energy. As we further increase the strength, one of these states becomes a true 
bound state after passing through zero at potential strength A^ = 1.55204. The other 
point that emerges from the bifurcation is a true virtual state that moves down the 
negative real axis as a function of A^. 

When the resonances are calculated by ECS (shown as red dots in figure [8]), 
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Figure 6. (color online) Continuation curves (A c = 0.5) shown in the full 
5R(n) X G(-!i) X A; space. Both start and end points of the curves are highlighted 
as points in respective colors. The four bifurcation points are indicated with big 
black dots. 



we do not resolve all these details. Especially the connections through the virtual 
states are missing from the picture. We discretize the two channel Hamiltonian on a 
finite difference grid with grid distance 0.03 and we use an exterior complex scaling 
transformation that starts at r — 12 and a rotation angle of ir/8. The exterior domain 
extends to 3?(r) = 15.6. 

As expected, at zero potential strength the eigenvalues of the exterior complex 
scaled Hamiltonian are the discrete eigenvalues of the kinetic energy operator. These 
are related to standing waves on the exterior complex scaled domain. These discrete 
continuum states will be rotated over twice the ECS rotation angle. Therefore, for 
each threshold we have a series of rotated eigenvalues. 

As the potential strength increases, these continuum states are attracted by the 
potential and become bound. The first bound state is formed from the smoothest 
continuum state associated with the first threshold. It becomes bound at A; « 0.5. 
This is in contrast to the numerical continuation results where this bound state 
originates from a virtual state. The difference between the numerical continuation 
and ECS is most prominent at small potential strengths. 

The second bound state in ECS is formed at potential strength A, ~ 1.57. 
Although this state should find its origin at zero potential strength in the smoothest 
continuum state of the second threshold, an avoided crossing appears instead. The 
curve originates from the second smoothest state of the first threshold and passes 
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Figure 7. (color online) The connectivity diagram of the different continuation 
curves. Curves only intersect at the four points bpi. Dashed (red/green) line 
connects two states from different channels through bp2 . See figure If 01 for a 
more detailed view of the connections. The end points of the curves are labeled 
accordingly but do not play a major role and their labels are omitted in other 
figures. 



Label 




Ai = A 2 


u 






E 


bpi 


2 


0852303e+00 


-1.4443524e+00 


-7 


0688100e- 


-02 


bp2 


1 


9571562e+00 


6.6636568e-01 


-8 


7009530e- 


-02 


bp3 


1 


5436785e+00 


8.8709701e-01 


-7 


2105286e- 


-03 


bp4 


4 


0009060e-03 


-8.6796523e-01 


-1 


0092983e- 


-02 



Table 2. Numerical values of the branch points of the resonant trajectories shown 
in figure [5] The coupling strength is A c = 0.5. 



through an avoided crossing with the trajectory of the state that starts from the 
smoothest state of the second threshold. The former curve partly represents the bound 
state, whereas the latter curve partly follows the trajectory of the Feshbach resonance 
between potential strengths 0.46 < < 1.45. However, as the curves are disconnected, 
this calculation misses the bifurcation at A,; = 1.54368 and the subsequent transition 
through a virtual state. Note that these issues occur near the line of ECS eigenvalues, 
starting from the first threshold going twice the ECS angle downwards the complex 
plane. Similary, the deviation of the trajectory around Ai = 0.46, as the Feshbach 
resonance turns into a virtual state, starts near an analogous line originating from the 
second threshold. These lines are drawn dashed in figure [5] 
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(b) fH{E) x 9f(E) x Aj space 



Figure 8. (color online) A comparison of the results of numerical continuation 
with the results of exterior complex scaling. The vertical axis shows the potential 
strength of the model problem while the real and imaginary parts of the energy 
are shown on the other axes. The blue and the green curves are the translation 
of the same curves in the u-plane from figures [5] and \E\ The red dots are the 
relevant ECS eigenvalues calculated for a range of A; values. The results of the 
two methods differ significantly in the regions where a resonance becomes a bound 
state. Also for small potential strengths we see significant deviations since ECS 
cannot resolve the virtual states. 
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Figure 9. (color online) SR(m) x projections of the continuation curves for 
different coupling strengths A c . Starting values are taken from table [T] 
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(c) A c = 0.3 (d) A c = 0.5 



Figure 10. (color online) Closeup view of the continuation curves of C2no, 
C2Uq and c\n\ obtained for different values of the coupling strength A c . As the 
coupling strength increases an interesting effect of reordering of the connections 
can be observed. Initially the trajectories do not interact. With slight coupling, 
two additional branch points appear connecting both curves. As the coupling 
increases, the branch points collide and disappear disconnecting the reordered 
curves. 

7. Conclusions and Outlook 

An automatic, robust and inherently efficient method for tracking parameter 
dependent resonant solutions of the Schrodinger equation is desirable since they play 
an important role in many quantum mechanical systems. Numerical continuation 
based on pseudo-arclength continuation has a proven track-record in being reliable 
and has been used in the study of various dynamical systems. 

We have shown earlier that numerical continuation can be applied to track bound 
and resonance states in a single and coupled channel Schrodinger equation with equal 
thresholds. However, when the thresholds arc different additional branch cuts appear. 
This leads to numerical difficulties. 

In this article we have shown that these numerical difficulties can be avoided when 
a suitable uniformization is applied. The channel momenta (i.e. the wave numbers) are 
then written as a function of a complex valued parameter. The numerical continuation 
is applied to this parameter combined with the variable parameters of the problem. 
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Unfortunately, this approach cannot be extended easily to systems with more than 
two channels because the uniformization procedure becomes too complex. 

We have applied the method to a two channel problem with Gaussian potentials 
and continued in the strength of the potential for various choices of coupling strength. 
Several branch points were detected and the continuation automatically identified the 
other branches emerging from these points. Transitions between bound and resonant 
states are easily taken by this method. In a similar way we could have continued in 
the threshold values or any other parameters starting from any solution point. 

The comparison of the results with an exterior complex scaling calculation shows 
significant differences in parameter ranges where the resonance transitions to a bound 
state. The numerical continuation results predict that this transition happens through 
a virtual state for the model problem. ECS, however, cannot resolve these virtual 
states and the resonance transitions directly to a bound state. 

In our calculations we have detected several branch points where different states 
meet. All the branch points we have identified, however, are bifurcations on the 
negative real energy axis where two virtual states meet. This occurs when a resonance 
becomes bound through a scenario that was already discussed by Nussenzveig [29] . 

In the future we will extend the method to higher dimensional problems with 
multiple reaction coordinates. 

Another possible future direction of research is to use two parameter continuation 
and automatically identify the exceptional points [30[ [31] [32] where two resonances 
coalesce at a complex valued energy that does not necessarely lies on the real axis. 
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